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\ The statistical mechanics of cold dark matter (CDM) particles interacting 

via a softened gravitational potential is reviewed in the microcanonical ensemble 
and mean- field limit. A phase diagram for the system is computed as a function of 
the total energy E and gravitational softening length e. For softened systems, two 
> \ stable phases exist: a collapsed phase, whose radial density profile p{r) is a central 

^ ' Dirac cusp, and an extended phase, for which p{r) has a central core and p(r) ~ 

^ ' ^-2.2 large r. It is shown that many A^-body simulations of CDM haloes in the 

Q ■ literature inadvertently sample the collapsed phase only, even though this phase 

is unstable when there is zero softening. Consequently, there is no immediate 
reason to expect agreement between simulated and observed profiles unless the 
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Qh' gravitational potential is appreciably softened in nature. 

O 

^ \ Subject headings: dark matter — galaxies: kinematics and dynamics — galaxies: 

^ ' structure — gravitation 

> ■ 

; 1. Introduction 

Cold dark matter (CDM) theory successfully describes many aspects of the formation of 
large-scale structure in the universe (Peebles 1982; Davis et al. 1985). However, mismatches 
do exist between its predictions and observations, such as the cusp-core controversy, missing 
satellites (Klypin et al. 1999; Moore et al. 1999a) and the angular momentum problem 
(Navarro & Benz 1991; Thacker & Couchman 2001). In particular, the cusp-core issue has 
provoked much debate. CDM simulations consistently yield density profiles with steeper 
inner slopes (power-law exponent between —1 and —1.5) than observational studies which 
have found a range of slopes, including constant density cores in dark matter dominated 
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low surface brightness galaxies (de Blok, McGaugh & Rubin 2001) and shallow slopes in 
clusters with gravitationally lensed arcs (Sand et al. 2004). These results, among others, have 
initiated discussion about the role baryons play in softening simulated cores (Athanassoula 
2004; Shen & Sellwood 2004) and the observational effects that may mask cusps in low surface 
brightness galaxies (de Blok, McGaugh & Rubin 2001; Swatcrs et al. 2003). Recent A^-body 
results demonstrate that, at the current resolution of simulations, the central power- law 
exponent does not converge to a universal value (Navarro et al. 2004). 

In numerical simulations, a softened gravitational potential is used to prevent the macro- 
particles (lO^-lO'^Afo) from experiencing artificially strong two-body interactions (Navarro, 
Frenk & White 1996, for example). The softening length, e, is chosen to maximise the res- 
olution while suppressing two-body effects over the simulation running time. In view of the 
ongoing disagreement regarding the form of p(r), it is important to clarify analytically, by 
a code- independent argument, whether the choice of e affects the physics of the system and 
hence p{r). In this paper, we employ the framework of statistical mechanics (Padmanabhan 
1990), drawing upon recent results on phase transitions in iV-body systems with attractive 
power-law potentials (Ispolatov & Cohen 2001; de Vega & Sanchez 2002). Self-gravitating 
particles behave qualitatively differently to many other statistical systems because gravity is 
an unscreened, long-range force. They are best examined within the microcanonical ensem- 
ble, where the energy and number of particles are fixed and phases with negative specific heat 
are allowed. We apply the results of the classical theory of self-gravitating, A^-body systems 
to demonstrate the effects of introducing a short distance cutoff in numerical simulations, in 
particular the effect on stability. 

The study of the thermal stability of self-gravitating systems has a long history. Lynden- 
Bell & Wood (1968) showed that spherical systems of point particles in a box with reflecting 
walls are gravitationally unstable below a critical temperature, collapsing catastrophically 
to a central point. Aronson & Hansen (1972) generalised this work to a spherical system of 
A^ classical hard spheres in contact with a heat bath, showing the gravothermal instability 
to be a general feature of self-gravitating systems held at a constant temperature. Hertel & 
Thirring (1971), investigating point fermions obeying the Pauli Exclusion Principle, showed 
that a stable low temperature phase can exist if the gravitational potential is softened, 
transforming the gravothermal instability to a phase transition to the low temperature phase. 
In this paper, we extend these results and use them to reinterpret some of the ambiguous 
results of numerical simulations of CDM haloes discussed above. A detailed comparison with 
preceding analytic work is presented in Section 3.2. 

Section 2 briefly reviews the formalism for treating A^ self-gravitating, coUisionless par- 
ticles statistically. In Section 3, we apply the formalism to compute p(r) analytically as a 
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function of E, the total energy, and e. The result is a thermodynamic phase diagram that 
contains both collapsed and extended haloes. In Section 4, we locate published A'"-body sim- 
ulations on the phase diagram and show they are biased towards the collapsed phase. This 
phase is unstable for e = 0, suggesting that collapsed haloes are an artificial by-product of 
the softened potential; there is no immediate reason to expect agreement between simulated 
and observed profiles unless the gravitational potential is appreciably softened in nature. We 
emphasise at the outset that it is not our intention to reproduce realistic CDM haloes with 
non-zero angular momentum and hierarchical clustering; rather, we demonstrate in a code- 
independent manner how the softening used in A'"-body simulations may artificially alter the 
density profiles found. 



The properties (e.g. energy, entropy) and collective behaviour (e.g. gravothermal catas- 
trophe) of a self-gravitating gas of CDM particles in thermodynamic equilibrium take differ- 
ent values when computed in different statistical ensembles because the long-range nature 
of the gravitational potential renders the system inseparable from its environment (Pad- 
manabhan 1990; Ispolatov & Cohen 2001; de Vega & Sanchez 2002). In this paper, we 
follow previous studies by considering the self-gravitating gas in the microcanonical ensem- 
ble (MCE), whose features are constant energy, volume and particle number. Particles do 
not evaporate from the system over time and the walls of the container are perfectly reflect- 
ing. The MCE is more appropriate than the canonical ensemble (CE) for three reasons: (i) 
it is unclear how to construct an external heat bath (required by the CE) for a long-range 
potential, because the system interferes with the environment (Huang 1987); (ii) states with 
negative specific heat are inaccessible in the CE (Padmanabhan 1990); and (iii) the equilib- 
rium density profile in the violently relaxed (Smoluchowski) limit is the singular isothermal 
sphere in the CE, contrary to observations (Sire & Chavanis 2002). 

The density of states, g{E). is the volume of the (6A^— l)-dimensional surface of con- 
stant energy E in phase space (xi, x^v, Pi, Pat), where (xj, pj) are the co-ordinates and 
momenta of the i-th particle. At any one moment, the system occupies one point in the 
6A^-dimensional phase space. For particles of equal mass m, one has 



2. Statistical mechanics of N self-gravitating pcirticles 



2.1. Density of states in the microcanonical ensemble 




(1) 
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where the first and second sums give the kinetic and potential energy, and the integral is 
over phase space volume. The gravitational potential, V, is given hy V = — Gm^lxj — Xj|~^ 
or, if the potential is artificially softened over a characteristic length (,hy V — — Gm^[(xj — 

The thermodynamic entropy S (up to a constant) and the temperature T of the system 
are defined in terms of g{E): 

S{E)^kBlng{E\ (2) 

^'-)^^^^- 

These quantities are hard to interpret when assigned to a system far from equilibrium. Note 
that g{E) diverges for C = and N > 2; any two particles can be brought arbitrarily close 
together, liberating an infinite amount of potential energy, so that the co-ordinate space 
integral diverges (Padmanabhan 1990). This is a serious problem because it is impossible 
to achieve thermodynamic equilibrium if g{E) diverges; the system does not have time to 
sample the infinite number of possible microstates with equal probability (Chabanol, Corson 
& Pomeau 2000). If the dark matter particles are fermions, the Pauli Exclusion Principle 
does prevent this problem. However, the fraction of the phase space volume sampled by 
N mildly relativistic CDM particles in a time t, given by {ct/ Rf^ {Cc^ /2GmN'^/^f^/^ , is 
exceedingly small for most proposed CDM particles, e.g. C^/m ~ 8xlO~(^^~^^W^/GeV for 
self- interacting dark matter (Spergel & Steinhardt 2000). 



2.2. Integral equations for p{r) in the mean-field limit 

The density of states is evaluated in the continuum (mean-field) limit by integrating 
over momentum and then expressing the remaining configurations as a functional integral 
over possible density profiles p{x) (de Vega & Sanchez 2002; Ispolatov & Cohen 2001), 

^ / C S C ^ 

where the effective dimensionless action 



.(A«.7.« = «+§//^g^^rf'x.d'x. (5) 
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^ I p(x)(i'^x — 7 — - ln/3 — y p(x) lnp(x)(i'^x 

and dimensionless energy 



C = ER/GM'^ 



(6) 



are quantities defined by Ispolatov & Cohen (2001). In (4) and (5), and throughout the 
remainder of this paper, the density profile and position co-ordinates x are written as di- 
mensionless quantities, relative to the total mass M — Nm and outer radius R of the system, 
with X — |x|/i? and e — (/R. 

Upon evaluating the functional integral by a saddle point method (which involves ex- 
tremising the action), (4) reduces to three coupled integral equations describing the density 
profile p{x), the central density po and the inverse temperature /3 = I/ZcbT. For a Newtonian 
potential, one has 



p{x) = Po exp 



1 

Po Jo 
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X 



iirxl exp 



/ p{xi)xi{\x + Xi\ — \x — Xi\)dxi 
Jo 

/ p{xi)xi{\x2 + Xi\ - \X2 - Xi\)dXi 
Jo 
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X2 



dX2, 



2/3 



^ + 47r^ / / p{xi)p{x2)xiX2{\xi + X2\ - \xi - X2\)dXidX2. 
Jo Jo 



(7) 
(8) 
(9) 



Note that the factor Anx^ in (8) was omitted due to a typographical error by Ispolatov & 
Cohen (2001). The solutions to these equations describe the density profile, entropy and 
temperature of an equilibrium system for a given energy. To obtain analogous equations for 
the softened potential, we replace |...| everywhere with [(...)^ + 6^]-*^/^. This extension is valid 
in the mean-field formalism for e small: correction terms are 0(e) (de Vega & Sanchez 2002). 

The mean-field limit is only meaningful physically for €7^0, otherwise g{E) diverges, de 
Vega & Sanchez (2002) verified the mean-field results against Monte-Carlo simulations and 
an alternative analytic method known as the Mayer cluster expansion, where the density 
of states is expanded as a combinatorial series in the dilute limit {C, ^ 1). The different 
approaches are in accord in the dilute limit. In order to verify the mean-field approach in the 
high-density hmit (^ <^ 1), relevant to the extended and collapsed phases studied here, we 
need to calculate the correction terms in this hmit, following de Vega & Sanchez (2002) — a 
project outside the scope of this paper. Nevertheless, to give a rough idea of these corrections, 
we note (by analogy) that they are of order 0{rfe, rje'^) in the CE, where 77 = Gm^N/RT in 
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the CE is a proxy for in the MCE. Note that the classical thermodynamic limit {N/V 
constant as N,V ^ oo) does not apply for gravitating systems; thermodynamic quantities 
are finite if proportional to N/V^/^ as N,V ^ oo. 

3. Radial density profile of a CDM halo 

We solve (7)-(9) for the radial density profile p{x) by the following iterative relaxation 
scheme (Ispolatov & Cohen 2001): given the current iterate of the profile, apply 

(9), (8) and (7) to compute f3'''^^^\ p^^^^ and p*{x) in that order, then apply p^*"*"-^^ (x) = 
ap*{x) + (1 — (t)p^'^\x) until the convergence criterion 

47r [ \p^'^^\x) - p^^{x)\x^dx < 5, (10) 
Jo 

is satisfied. We typically adopt 6 = lO^*' and 0.01 < a < 1 (5 ^ o") in this work. The 
softening can be introduced into this scheme in two ways: (i) as a nonzero lower limit of 
integration in the integrals in (7)~(9); and (ii) in the potential, V = — Gm^[(xi — x^)^ + 
^2j-i/2 gQ^]^ approaches were tested and found to produce quahtatively similar behaviour; 
we concentrate on the latter in this work as it is more closely allied to A'"-body simulations. 

3.1. Stable versus unstable phases: e = 

With no softening present in the gravitational potential, a stable solution of (7)-(9) 
formally exists above a cutoff energy ,^ > ~ —0.335. The density profile of the halo exhibits 
a flat central core, with dp/dx — > as x — > 0, and near-isothermal wings, with p{x) oc 
(a ~2.2) as X — > oo, as illustrated in Figure 1(a). We refer to it as the extended phase. It 
agrees with the solution for secondary infall onto a spherical perturbation (Bertschinger 1985) 
and behaves asymptotically like the spherical, thermally conducting polytrope (Lynden-Bell 
& Eggleton 1980) and infinite-dimensional Brownian gas (Sire & Chavanis 2002). 

For ^ < and e = 0, a formal solution of (7)-(9) does not exist. The entropy and 
temperature are discontinuous at this cutoff energy as shown in Figure 1(b). This is the 
well-known gravothcrmal catastrophe (Antonov 1962). Note that, for ^ = —1/4, the singular 
isothermal sphere p{x) = (47rx^)^^, po = (47re^)~^ and /9 = 2 is always a solution of (7)- 
(9), as can be verified analytically, but it is not stable and so the iterative procedure never 
converges to it, but rather to Figure 1(a). 
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3.2. Collapsed versus extended phases: e 7^ 

If the gravitational potential is softened, a stable phase exists for all values of ^. For 
^ > ^c'^^ ~ 0, the halo is extended as in Figure 1(a), with a flat core and near-isothermal 
envelope, p{x) oc x^"^-^. However, for ^ < {c, the halo is collapsed. Figure 2(a) displays the 
density profile of such a collapsed halo for e = 10~^ (note that both axes are logarithmic). 
The halo has a steep Dirac peak ('cusp') at x — 0: p{x) is flat for x < e, decreases as a large 
inverse power of x for x > e and flattens for x — > 1. For intermediate energies in the range 
< ^ < ^c^\ the system is bistable: the halo can be either extended or collapsed depending 
on the initial conditions and the route to equilibrium. Figure 2(b), a plot of entropy and 
temperature as a function of energy, illustrates this bistability and the hysteresis to which it 
can lead. S{^) and /3{^) jump discontinuously at both .^c''''' and ^c- If ^ enters the intermediate 
range from below, the halo remains collapsed until ^ exceeds ^c^'*. Alternatively, if ^ enters 
from above, the halo remains extended until ^ is reduced below (Ispolatov & Cohen 2001). 

The critical energy and the collapsed and extended profiles we obtain, are consistent 
with previous analyses (Antonov 1962; Aronson & Hansen 1972; Padmanabhan 1990). For 
example, Aronson & Hansen (1972) find a phase transition at ^ ~ —0.3, consistent with 
our value — —0.335, and a critical reciprocal temperature at the transition in the range 
Pc — 1.2 — 1.6 for e = lO"^'^, consistent with ~ 2 in this paper (e = 10~^). A more precise 
comparison is prohibited by the adoption of the CE rather than the MCE in most previous 
work. The van der Waals model proposed by Padmanabhan (1990) is an exception; it is 
examined in detail in the following section. The phenomenon of bistability was overlooked 
until the work of Ispolatov & Cohen (2001). 

The behaviour of the system depends somewhat on the choice of the relaxation param- 
eter (7, defined at the start of this section. Table 1 displays the minimum softening length 
for which the system makes the transition to a stable collapsed phase, for a given value of 
a. The phase transition from the extended to the collapsed phase is increasingly delicate as 
e decreases; for larger values of a, the system is less likely to reach the critical point where 
the phase transition occurs and thus remains in the extended phase. Although the effect is 
numerical, it potentially reflects the relative hkelihoods of the possible routes that the real 
system can take to equilibrium. 

3.3. van der Waals equation of state 

An alternative, phenomenological way to model a CDM halo interacting via a softened 
gravitational potential is to solve the Lane-Emden equation for a nonideal, isothermal gas 



Fig. 1. — (a) Equilibrium radial density profile for a system with ^ > a-nd e = 0. The 
profile is similar to a softened isothermal sphere, but with p{x) oc x~^^^ as x ^ 1. (b) 
Dimensionless entropy (monotonically increasing) and inverse temperature (peaked curve) 
as a function of energy for ^ > when the gravitational potential is not softened. Below 
^c— — 0.335, no stable solution exists. 




Fig. 2. — (a) Radial density profile for the collapsed phase of a potential with softening 
e = lO"*^. (b) Dimensionless entropy (thick curve) and inverse temperature (thin) versus 
energy. The collapsed and extended phases are labelled C and E respectively. Both quantities 
display discontinuous jumps at the phase transitions at — —0.335 and ^c^^ ~ 0.0. If the 
energy is increased from below ^c, the system remains in the collapsed phase until ^ = ^c^\ 
when it jumps to the extended phase. If the energy is reduced from above the system 
remains in the extended phase until ^ = ^c, when it jumps to the collapsed phase. 
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obeying a van der Waals equation of state, P oc pT{l — p/ Pm)~^, where P denotes the 
pressure and Pm is the maximum density allowed by hard- sphere packing of the CDM particles 
(Aronson & Hansen 1972; Padmanabhan 1990). The analogy with a softened gravitational 
potential implies Pm ~ me~^, although it is clear from the outset that this analogy is inexact; 
the hard-core, van der Waals potential effectively excludes a volume ~ around each 
particle, whereas softened gravity suppresses the mutual acceleration of particles separated 
by a distance e without preventing them from 'coasting' even closer together. We compare 
our solutions of (7)-(9) with the van der Waals model by integrating the nonideal Lane- 
Emden equation 

1 d 



^2 flj. 

from r — Ri to r — R, obtaining 



p dr V 1 - p/Pn 



1 dp 
p{R) 'dr 



( > = -4.Gp (11) 



= -I3M^/M , (12) 



^1 dp 

^~ i?V(^l)[l-p(i?l)/pJ2 ^ 



where Mi denotes the mass enclosed in the volume Ri < r < R. [We use p{R) <^ pm to 
simplify (12).] The quantities p(-R), p(-Ri), their derivatives. Mi, and (3 can be extracted from 
our numerically computed profiles, which satisfy (7)-(9), in order to find pm as a function of 
^ and e and hence compare with the results of Aronson & Hansen (1972) and Padmanabhan 
(1990). The inner integration limit Ri must be positioned carefully near the inflection point 
of the central Dirac peak, e.g. at r = 10~^'^ in Figure 2, in order to avoid the innermost grid 
cells, where the numerical solution is noisiest, while ensuring that p{R\)/ p^ is not too small, 
to avoid roundoff error when solving (12) for p^. 

In Figure 3, we compare the density of states and van der Waals models by plotting 
l/j3 versus ^ for both models in the MCE. The open triangles indicate solutions of (11) 
for 10~^ < a = M/Airp^-aR^ < 10~^. The boxes and asterisks indicate solutions of (7)-(9) 
for e = 10~^ and 10~^ respectively, with a tolerance of 6 = 10~^. (We have verified that 
the results are unchanged for 10~^ < S < 10~^.) Applying the procedure in the previous 
paragraph to compute pm, we find logo = —10.1 and —9.00 for the boxes and asterisks 



^ f min 

OOl 10^ 

0.03 10"^' 

0.05 10-^ 

0.1 10-2 



Table 1: Minimum softening, emin, for which a stable collapsed solution is found for a given 
relaxation parameter a. 
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respectively. For ^ > ^c, the two models are in accord, as expected; the halo is extended, so 
the softening (and the precise value of a) are not important in the dilute regime p <^ pm- 
For ^ < ^c, the two models differ appreciably. The density of states calculation predicts less 
variation of T with E than the van der Waals model, for a given value of a. We confirm 
the trend, apparent in Figure 4.11 of Padmanabhan (1990), that T increases with \a\ when 
^ < is fixed. However, the trend is confirmed in the range — 10 < log a < —9, which does 
not overlap with the range —4 < logo < —2 investigated by Padmanabhan (1990). The 
effective value of pm predicted by the density of states calculation is systematically greater 
than anticipated in van der Waals models published previously. 

3.4. Phase diagram on ^ — e plane 

We can produce a phase diagram for the system by plotting the power-law exponent of 
the density profile as a function of ^ and e. The logarithmic slope of the density profile is 
defined as 

, , c?ln[p(a;)l , , 

We evaluate the logarithmic slope in the inner halo, at x = e, and also at the box edge, 
X = 1, to consistently characterise the two phases. Figures 4(a) and 4(b) display the phase 
diagrams obtained from p{e) and p{l) for a = 0.01; note that varying the value of a does not 
affect the results noticeably. The phase transition occurs where the contours are bunched, at 
^ ~ —0.335. For e < 10~^, p(e) is large (Dirac cusp) and p{l) is almost zero (flat envelope) 
in the collapsed phase. For e > 10~^, the system does not undergo a phase transition. The 
entropy and temperature are continuous for all ^ and the profile is extended everywhere. 

Interestingly, the central Dirac peak steepens as e decreases, in the collapsed phase. 
Figure 5 plots p{x) as a function of logarithmic radius for 10~^ < e < 10~^ and at a fixed 
energy ^ = —0.5. The location and amplitude of the maximum of p{x) is a strong function 
of e, with p{x) — > as X — > 0, 1 for collapsed haloes. The dynamic range of numerical studies 
generally extends from a few times e to the virial radius. In this range, p{x) is strongly 
affected by the softening length used. The classical NFW profile (Navarro, Frenk & White 
1996) and the profile of the extended phase at > .^c are also plotted for comparison. Clearly 
the NFW profile, although it is cuspy, does not match exactly the collapsed (or indeed the 
extended) phase profile. Wc are unable to explain unambiguously why the NFW profile is 
not reproduced, but remind the reader that our analysis neglects several effects, such as 
hierarchical clustering, non-zero angular momentum, and cosmological expansion, which are 



- 11 - 



o 
O 



2.0 
1 .5 

1 .0 

0.5 

0.0 

0.5 
1 .0 



A 



A 




12 3 

ER/GM' 



Fig. 3. — Temperature T as a function of total mechanical energy E for a range of softening 
lengths e. The open squares and asterisks denote CDM haloes with e = 10~^ and 10"^ 
respectively, modeled by the density of states formalism (7)-(9). These haloes correspond to 
log a = —10.1 and —9.00 in the van der Waals formalism. The open triangles denote CDM 
haloes with —4 < log a < —2, modeled by the van der Waals formalism (Padmanabhan 
1990). For any given E in the plotted range, T is the same for all values —4 < log a < —2; 
moreover, these a values are systematically greater than the predictions of the density of 
states formalism. 
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present in full A^-body simulations. The latter effect is discussed in more detail below. 

3.5. Cosmological expansion 

Our calculations are performed in a static background spacetime. In the central parts 
of the halo, where we focus most of our attention, CDM particles are tightly bound and 
effectively decoupled from the cosmological expansion. Consequently, the phase transition 
at is essentially unaffected, as long as the central volume where the cosmological expansion 

can be neglected is still large enough to encompass the 'core-halo' structure (centre-to-edge 
density contrast > 709) required for the gravothermal catastrophe to occur (Lynden-Bell & 
Wood 1968). From this perspective, a static background metric is a good approximation. 

In the outer halo, CDM particles are loosely bound and p{x) may be determined partly 
by the expansion and certain specific cosmological parameters. For example, there exists 
a mapping between the linear and evolved two-point correlation functions of haloes which 
implies power-law density profiles with exponent — 3(n + 4)/(n + 5), where n is the index of 
the power spectrum of initial fluctuations (Hamilton et al. 1991; Padmanabhan & Engineer 
1998). By the same token, simulations suggest that the dependence of the shape of the 
profile on initial cosmological parameters is weak; this empirical result holds for SCDM 
(Qm = 1-0), LCDM (Qm — 0.3, = 0.7) and open {Qm — 0.3) cosmogonies, as well as a 
wide mixture of hot and cold dark matter components (Huss, Jain & Steinmetz 1999a). Note 
that the origin of the NFW scale radius observed in simulated haloes, r^, is not explained 
by the (self-similar) Padmanabhan & Engineer (1998) mapping. However, the total power 
spectrum, P{k), is proportional to the product of the Fourier transformed halo profile and the 
power spectrum of the distribution of halo centres, Pccnt(^) (Padmanabhan 2002), implying 
p{x) oc at large k [deep minima of the gravitational potential, V, with Pcent(^) oc Pv{k)], 
and p{x) (X x~^ at small k [quasi-linear regime, Pccnt(^'") P{k)] (Padmanabhan 2002). 

Huss, Jain & Steinmetz (1999b) demonstrated empirically that cosmological expansion 
does influence p{x), but that it is one of several relevant factors. They simulated a halo 
with all tangential components of the gravitational force artificially set to zero, reproducing 
the spherical infall solution with a power- law exponent of ~ —2.2 and no break in slope 
(Bertschinger 1985). They also showed that is not determined solely by expansion, found 
evidence for the importance of angular momentum in the system, and concluded that the 
unbroken power-law profiles of van Albada (1982) in the absence of expansion were inade- 
quately resolved for the purpose of testing the existence of r^. Of course, the scale radius 
(and hence the concentration parameter) of a halo does depend critically on its formation 
epoch: is fixed by the overdensity, Sc, which is, in turn, a function of the collapse redshift 



(Navarro, Prenk & White 1997). 



4. Compctrison with numerical simulations of CDM haloes 

The key result from Section 3 from the perspective of A'"-body simulations is that any 
simulation with e 7^ can produce stable haloes at energies ^ < that do not yield stable 
haloes in true (e = 0) gravity. Furthermore, these stable e 7^ haloes are collapsed, whereas 
stable e = haloes are extended. In this section, we show that many pubhshed A'"-body 
results inadvertently sample the collapsed phase only. We confine ourselves to studies that 
report the specific values of ^ and e explored. 



4.1. Softening length 

Simulations have been performed over a range of N and with a range of resolutions in an 
attempt to place bounds on the optimum softening length (Ghigna et al. 2000; Splinter et al. 
1998; Moore et al. 1998). van Kampen (2000) argued that the choice e ^ 0.5ri/2N~^^^, where 
ri/2 is the half-mass radius of the system, strikes a balance between too short a relaxation 
time (e too small) and excessive particle clustering (e too large). Similarly, Athanassoula 
et al. (2000) find an optimal length e = 0.32A^~°'^^ for a 7 = Dchncn sphere. These 
different criteria define a range of e within which most modern simulations are performed. 
The softening length in high resolution simulations is 1-5 kpc for galaxy-sized haloes (Reed 
et al. 2003), corresponding to e ~ 10~^ if it! is taken to be the virial radius. The studies 
investigated below occupy the range 10~^ < e < 10~^. 



4.2. Total energy of a halo 

The total energy of a simulated halo is rarely quoted in published studies. We therefore 
calculate ^ from the quoted mass, size and concentration parameter, c = rs/r2oo, where 
is the characteristic radius of the halo and r2oo = -R is the radius at which the halo density 
has dropped to 200 times the background. For a classical NFW halo, the kinetic energy, K, 
is given by (Mo, Mao & White 1998) 



Knfw = (14) 
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f{c) 



c[l - 1/(1 + c)2 - 21n(l + c)/(l + c)] 



2[c/(l + c)-ln(l + c)]2 
assuming circular orbits. The potential energy is (Binney & Tremaine 1987, eq. 2P-1) 



(15) 



WFW — 



2R 



-h(c), 



h{c) 



[ x~'^[{l + cx)~^ -l + ln{l + cx)]^dx, 
Jo 



(16) 
(17) 



assuming spherical symmetry. Similarly, for the halo profile found by Moore et al. (1999b), 
we obtain 



Km 
Um 
9{c) 



AR 

GM^ 



1 + 



9ic) 



In^ (1 + c3/2) 
9{c) 



2R ln'(l + c3/2)' 

Jo 



(18) 
(19) 
(20) 



Wc estimate the total energy, E, from (14)-(20) in three ways: (a) E^"^ = -K, which 
assumes the virial theorem and circular orbits (as stated above); (b) E^''^ = U/2, which 
assumes the virial theorem only; and (c) E^'^^ — K + U, which assumes circular orbits, but 
not the virial theorem. All three approaches take p(x) to be spherically symmetric. 



4.3. Position on the phase diagram 

Wc now locate on the phase diagram some examples of published A^-body simulations of 
dark matter haloes in a ACDM cosmology, for a range of halo masses from dwarf galaxy (~ 
IO^^Mq) to cluster (~ IO^^Mq) size. Prom the output parameters, including c, we calculate 
^ and e (normalised by r2oo)- Table 2 summarises these data and the energy estimates 
obtained by the three methods (a)-(c) in Section 4.2. In the final column, comments are 
made identifying the specific haloes chosen from the referenced work. In Figure 6, we place 
the simulated haloes on the ^'^''-'-e phase diagram [Figure 4(a)]. ^^''^ is preferred as it does 
not rely on the restrictive assumption of circular orbits. 

It is striking that the pubhshed haloes exist exclusively in the collapsed phase, near the 
left-hand edge of the diagram, or else in the bistable region, where stable solutions exist for 
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both collapsed and extended phases. They do not exist in the extended phase at the right- 
hand edge of the diagram. Any NFW halo with energy ^^"'^ necessarily lies in the collapsed 
regime, because one has /(c) > 2/3. Alternatively, an NFW halo with energy ^^'^^ lies in the 
collapsed regime unless one has c < 8 (corresponding to more massive systems). Similarly, 
all Moore haloes have energy ^^""^ < —3/8, and are restricted to the collapsed regime for all 
c > 0, while a Moore halo with energy ^^''^ is collapsed unless one has c < 4. In other words, 
if the total energy is estimated from the kinetic energy and virial theorem, both the NFW 
and Moore profiles are in the collapsed phase. Otherwise, small concentration parameters 
allow either phase (in the bistable region), depending on the detailed route to equilibrium. 
Without the virial assumption, haloes with energy ^^'^^ lie in the collapsed or bistable regions 
as well, mostly in the former. 



5. Conclusion 

We have investigated the equihbrium configurations of N self-gravitating coUisionless 
particles, interacting via a softened gravitational potential, in the MCE and mean-field limit. 
Below a critical energy, — —0.335, a system with e ^ exists in a stable, collapsed phase. 
This phase is unstable for pure gravity (e = 0). Above another critical energy ^c'^^ ^ 0, both 
softened and unsoftened systems exist in an stable, extended phase. In the intermediate 
region .^^ < ,^ < ^c^\ both the collapsed and extended phases are accessible; the detailed 
route to equilibrium determines which one is picked out. The density profiles for the extended 
and collapsed phases are qualitatively different. The extended profile has a flat core and near- 
isothermal outer envelope. The collapsed profile is a centrally condensed Dirac peak, whose 
logarithmic slope depends on e. 

We compare our results with published A^-body simulations by using the softening 
parameter, e, size, r2oo, and concentration parameter, c, to place simulated haloes on the ^-e 
phase diagram. We find that many published simulations inadvertently sample the collapsed 
phase only, even though this phase is unstable for pure gravity and arguably irrelevant 
astrophysically. 

We remind the reader that we neglect several effects that are important in real CDM 
haloes, such as hierarchical clustering, nonzero angular momentum, and cosmological ex- 
pansion. Our results elucidate some of the artificial behaviour that a softened potential can 
introduce; they are not a substitute for a full A"-body calculation. 

We are grateful to Bruce McKellar for extensive discussions on the theoretical basis of 
the mean-field equations, and the anonymous referee for useful comments that improved the 
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Fig. 4. — Contour plot of the logarithmic slope , \p\, of the density profile in its (a) inner, 
X = e, and (b) outer, x = 1, regions, as a function of energy, ^. The phase transition is clear 
at ^ ~ —0.335 for e < 10~^. The inner slope \p{e)\ steepens as e is reduced. For e > 10^^, 
the system does not collapse. 



Study 


M (10^2^0) 


r2oo (kpc) 


c 




^(^) 




e 


Comment 


Ref. 


NFW 


2.9 


172 


17.54 


-0.76 


-0.52 


-0.28 


0.01 




1 




22.7 


733 


15.38 


-0.73 


-0.48 


-0.23 


0.01 




1 


Hubs 


500 


1360 


6.3 


-0.55 


-0.29 


-0.03 


0.0037 


CDM run 


2 


Moore 


430 


1950 


4 


-0.57 


-0.32 


-0.07 


0.0005 


cluster 


3 


Reed 


0.188 


119 


28 


-0.93 


-0.68 


-0.43 


0.0020 


dwfl 


4 




40 


705 


12.5 


-0.68 


-0.43 


-0.18 


0.0009 


grpl 


4 


Hayashi 


2.2 


212.7 


5.3 


-0.52 


-0.27 


-0.02 


0.0021 


G3/2563 


5 



Table 2: Recent A'"-body simulations of CDM haloes. The mass, M, size, r2oo, concentration 
parameter, c, and softening length, e, are measured from the output. Equations (14)-(20) 
are used to calculate the halo energy in three ways, and ^('^l All calculations assume 

an NFW profile, except for Moore et al. (1999b) which uses the Moore profile. 



References. — (1) Navarro, Frenk & White (1996); (2) Huss, Jain & Steinmetz (1999b); (3) Moore et al. 
(1999b); (4) Reed et al. (2003); (5) Hayashi et al. (2003) 
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Fig. 5. — Logarithmic slope, p, as a function of logarithmic radius for ^ = —0.5 and e = 10~^ 
(long dash), 10~^ (dash-dot-dot-dot), 10"^ (dash-dot) and 10^^ (dotted). The slope is a 
maximum at x slightly above e. Also plotted is an NFW profile (solid) and an extended 
phase profile (short dash, ^ = —0.3). 







— 2.5 




-0.2 



Fig. 6. — Energy. ^^''^ . and softening length, e, of published A^-body CDM haloes overlaid on 
a contour plot of j:i(e), as summarised in Table 2. The asterisks, plus signs, triangle, box and 
diamond denote the results of Navarro, Frenk & White (1996), Reed et al. (2003), Hayashi 
et al. (2003), Moore et al. (1999b) and Huss, Jain & Steinmetz (1999b) respectively. 
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